#include<stdio.h>
#include<gmp.h>
#include<math.h>
#include<pthread.h>

#define nIterations 13 			//constante que define o numero maximo de iterações a fim de se obter 10 milhoes de digitos

pthread_t threads[2];			//declaracao do numero maximo de threads
pthread_mutex_t mutex[nIterations];

mpf_t y_k[nIterations], y_k1; 		//declaracao de y_k e y_k1
mpf_t a_k[nIterations], a_k1;		//declaracao de a_k e a_k1

void *calc_y(void *arg){
	mpf_t tmp_1, tmp_2, tmp_3, tmp_4, tmp_sup, tmp_inf; 		//declaracao das variaveis temporarias
	mpf_inits(tmp_1, tmp_2, tmp_3, tmp_4, tmp_sup, tmp_inf, NULL); 	//inicializacao das variaveis temporarias
	
	int i;
	for(i = 1; i < nIterations; i++){
		/* Calcula o y_k+1.*/
		mpf_pow_ui(tmp_1, y_k[i-1], 4);
		mpf_ui_sub(tmp_2, 1, tmp_1);
		mpf_sqrt(tmp_3, tmp_2);
		mpf_sqrt(tmp_4, tmp_3);
		mpf_ui_sub(tmp_sup, 1, tmp_4);
		mpf_add_ui(tmp_inf, tmp_4, 1);
		mpf_div(y_k1, tmp_sup, tmp_inf);
		/* Fim. */
		
		/* Liberamos o uso de y_k. */
		mpf_set(y_k[i], y_k1); 	//armazena o valor de y_k1 em y_k
		pthread_mutex_unlock(&mutex[i]); //destrava o mutex, sinalizando o fim do calculo de y
	}

	mpf_clears(tmp_1, tmp_2, tmp_3, tmp_4, tmp_sup, tmp_inf, NULL); //limpa as variaveis temporarias
}

void *calc_a(void *arg){
	mpf_t tmp_1, tmp_2, tmp_3, tmp_4, tmp_sup, tmp_inf; //definicao das variaveis temporarias
	mpf_inits(tmp_1, tmp_2, tmp_3, tmp_4, tmp_sup, tmp_inf, NULL); //inicializacao das variaveis temporarias

	int i;
	for(i = 1; i < nIterations; i++){
		pthread_mutex_lock(&mutex[i]); //a thread espera o termino do calculo do y

		/* Calcula a_k+1. */
		mpf_add_ui(tmp_1, y_k[i], 1);
		mpf_pow_ui(tmp_2, tmp_1, 4);
		mpf_mul(tmp_sup, tmp_2, a_k[i-1]);
	
		mpf_pow_ui(tmp_2, y_k[i], 2);
		mpf_add(tmp_3, tmp_1, tmp_2);
		mpf_mul(tmp_4, y_k[i], tmp_3);
		mpf_set_ui(tmp_1, 2);
		mpf_pow_ui(tmp_2, tmp_1, 2*(i-1)+3);
		mpf_mul(tmp_inf, tmp_2, tmp_4);
	
		mpf_sub(a_k1, tmp_sup, tmp_inf);
		/* Fim. */

		mpf_set(a_k[i], a_k1); //armazena o valor de a_k1 em a_k[i]
	}
	mpf_clears(tmp_1, tmp_2, tmp_3, tmp_4, tmp_sup, tmp_inf, NULL); //limpa as variaveis temporarias
}	

int main(){
	int i;
	mpf_set_default_prec(35000000); //precisao padrao: 35 milhoes de bits (aprox. 10 milhoes em decimal)
	mpf_t pi;

	FILE *file; //arquivo de saida		

	file = fopen("PIP", "w+");              //abre arquivo PIP para escrita
	if(file == NULL){
		printf("Erro ao abrir o arquivo\n");
		return 0;
	}
	
	mpf_inits(pi, y_k1, a_k1, NULL); 	//inicializacao das variaveis 
	for(i = 0; i < nIterations; i++){
		pthread_mutex_init(&mutex[i], NULL);
		pthread_mutex_lock(&mutex[i]);
		mpf_inits(a_k[i], y_k[i], NULL);
	}
	pthread_mutex_unlock(&mutex[0]);

        mpf_sqrt_ui(a_k[0], 2); 	//
        mpf_mul_ui(a_k[0], a_k[0], 4);	// calculo do a0 = 6 - [(4)*(2^1/2)]
        mpf_ui_sub(a_k[0], 6, a_k[0]);	//

        mpf_sqrt_ui(y_k[0], 2);		// calculo do y0 = 2^1/2 - 1
        mpf_sub_ui(y_k[0], y_k[0], 1);	// 

	
	pthread_create(&threads[0], NULL, calc_y, NULL); //cria a thread que calcula y_k
	pthread_create(&threads[1], NULL, calc_a, NULL); //cria a thread que calcula a_k
	pthread_join(threads[1], NULL);			 //espera as threads executarem
	
	mpf_ui_div(pi, 1, a_k[nIterations - 1]); 	 //calcula o inverso do valor de ak1 encontrado, que sera o valor do PI
//	gmp_printf("Pi %.52Ff\n", pi); 			 //imprime o valor do PI na tela, com 52 casas decimais.
	mpf_out_str(file, 10, 0, pi); 			 //imprime o valor do PI no arquivo, com 10000000 casas decimais.	

	//fecha arquivo e limpa variaveis//
	fclose(file);
	file = NULL;	
	mpf_clears(pi, y_k1, a_k1, NULL);
	for(i = 0; i < nIterations; i++){
		pthread_mutex_destroy(&mutex[i]);
		mpf_clears(a_k[i], y_k[i], NULL);
	}
	//fim//
	
	return 0;
}
